The apriori estimate of the Lieb-Robinson lightcone

The approach presented above is useless from a practical point of view: we first solve the many-body problem, then estimate the lightcone. To make this approach useful, we need to reverse the steps: we somehow make an a priori estimate of the lightcone, then solve the many-body problem inside this lightcone.

Maximum entropy estimate of the Lieb-Robinson lightcone

The first approach to estimate the lightcone is to employ maximum entropy-like considerations. Let us repeat the illustrative Hamiltonian of atom inside a waveguide:

\[ H\left(t\right) = E_{at} \sigma_{z} + \sigma_{x} \cdot f(t) + h \sigma_{+} a_{0} + h a^{\dagger}_{0} \sigma_{-} + \sum\limits_{i=0}^{\infty} (\varepsilon a^{\dagger}_{i} a_{i} + h a^{\dagger}_{i+1} a_{i} + h a^{\dagger}_{i} a_{i+1}) \]

We see that the atom speaks to the waveguide by placing a quantum at the site 0. Therefore, in order to estimate the lightcone, let us see how a single quantum freely propagates from the site 0.

Here we have a one-particle first-quantized problem. The wavefunction \(\left| \phi\left(t\right) \right\rangle\) of a single quantum is a superposition of its position on different sites of the chain:

\[ \left|\phi\left(t\right)\right\rangle = \sum_{k=0}^{\infty} \phi_k\left(t\right) \left| k \right\rangle, \]

where \(\left| k \right\rangle\) represents the quantum which is localized on the \(k\)th site.

The initial value is

\[ \left|\phi\left(0\right)\right\rangle = \left| 0 \right\rangle, \]

or

\[ \phi_k\left(0\right) = \delta_{k0} \]

The free propagation of the quantum is given by a first-quantized tridiagonal Hamiltonian

\[\begin{split} H = \left[\begin{array}{cccc} \varepsilon & h & 0 & \dots\\ h & \varepsilon & h & \dots\\ 0 & h & \varepsilon & \dots\\ \vdots & & \vdots & \ddots \end{array}\right] \end{split}\]

Let us solve numerically for \(\left| \phi\left(t\right) \right\rangle\):

import numpy as np
import tools

tmax = 800
dt = 0.01
t = np.arange(start = 0, stop = tmax, step = dt)
nt = t.size

ns = 100 # number of chain sites we keep

h = 0.05 # chain hopping
e = 1.0  # chain on-site energy

H = tools.tridiag([e]*ns, [h]*(ns-1)) # Hamiltonian for free propagation of a single quantum

phi_ini = np.zeros(ns, dtype = np.cdouble)
phi_ini[0] = 1 # initially quantum is on the site 0

psi_lc = np.zeros((ns, nt), dtype = np.cdouble) # Here we store the propagated orbitals

def Ht(t):
    return H

for i, psi in tools.evolution(start_index = 0, end_index = nt, H = Ht, dt = dt, initial_state = phi_ini):
    psi_lc[:, i] = np.copy(psi)

We see that the single-quantum problem is easily solved for large times and large chain sizes.

Let us animate the motion of the wavepacket:

import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.animation import FuncAnimation

psi_grid = psi_lc[:, ::100]
x = range(ns)
t_ = t[::100]

def update_fi(frame, line, x):
    line.set_ydata( np.abs(psi_grid[:, frame])**2)
    return [line]
    
fig, ax = plt.subplots()
ax.set_xlabel("chain sites")
ax.set_ylabel("|psi(x)|**2");
y = np.abs(psi_grid[:, 0])**2
    
line, = ax.plot(x, y)

nt_ = len(t_)
    
time = np.arange(0, nt_, 1)

animation = FuncAnimation(
    fig,                
    func=update_fi,    
    frames=time,       
    fargs=(line, x),    
    interval=30,       
    blit=True,         
    repeat=True)       

from IPython.display import HTML
display(HTML(animation.to_jshtml()))
plt.close()

Once Loop Reflect

We see that the wavepacket of the quantum just flies away from the point of emission. We can intuitively illustrate this process as a ray originating at the space-time point \((x=0, t=0)\):

Here the ray corresponds to the whole wavepacket trajectory \(\left\{\left|\phi\left(t\right)\right\rangle: t \geq 0 \right\}\). It is created on site \(x=0\) via \(\widehat{a}_0^{\dagger}\). Then the wavepacket flies away according to the animation.

Now suppose that the atom emits quanta into the waveguide during some time interval \(\left[0, t \right]\). Then at each time moment new quantum is placed to the site 0, which becomes a source of a new ray. This leads to the following picture:

We see that by the time \(t\), there are just newly created quanta in the state \(\left|\phi\left(0\right)\right\rangle\). There are also the quanta which where emitted at time \(0\), so that by \(t\) they are in the state \(\left|\phi\left(t\right)\right\rangle\). There are also quanta which were emitted at the intermediate times, so that the whole range of states \(\left\{\left|\phi\left(\tau\right)\right\rangle: 0 \leq \tau \leq t \right\}\) appears at \(t\).

Then the maximum entropy method is that we do not want to consider what was the actual quantum state of the atom when these quanta were emitted. We also do not want to take into account the correlations between the quanta. We just want to take into acccount the fact that under the free propagation their possible single-particle states are covered by \(\left\{\left|\phi\left(\tau\right)\right\rangle: 0 \leq \tau \leq t \right\}\).

Each \(\left|\phi\left(\tau\right)\right\rangle\) corresponds to a pure density matrix \(\left|\phi\left(\tau\right)\right\rangle\left\langle \phi\left(\tau\right)\right|\). Then the maximum entropy estimate of the one-body reduced density matrix of all the quanta emitted during \(\left[0, t\right]\) is

(6)\[ \widehat{\rho}_{+}\left(t\right) =\mathcal{N}^{-1}\int_0^t d\tau \left|\phi\left(\tau\right)\right\rangle\left\langle \phi\left(\tau\right)\right|, \]

with some irrelevant normalization constant \(\mathcal{N}\).

Lieb-Robinson metric

Given the explicit expression (6), we can pose a number of interesting questions.

Suppose we have some arbitrary one-quantum state \(\chi\),

\[ \left|\chi\right\rangle = \sum_{k=0}^\infty \chi_k \left| k \right\rangle. \]

We may ask the question, whether \(\chi\) is statistically significant for the time evolution during \(\left[0, t\right]\). The average significance is then given by \(\left\langle \chi \right| \widehat{\rho}_{+}\left(t\right) \left| \chi \right\rangle\). We need a treshold in order to decide whether this state can be discarded. There is a one-quantum state \(\left|\phi_m\right\rangle\) which is the most significant on \(\left[0, t\right]\):

\[ \phi_m = \arg\max_{\phi} = \left\langle \phi \right| \widehat{\rho}_{+}\left(t\right) \left| \phi \right\rangle. \]

It’s average significance is given by the maximal eigenvalue \(\pi_m\left(t\right)\) of \(\widehat{\rho}_{+}\left(t\right)\) due to the Rayleigh-Ritz variational principle.

Then we choose a relative significance treshold \(r_{cut}\), e.g. \(r_{cut}=10^{-4}\) and introduce the Lieb-Robinson metric

\[ g\left(\chi, t\right) = \left\langle \chi \right| \widehat{\rho}_{+}\left(t\right) \left| \chi \right\rangle - r_{cut} \pi_m\left(t\right). \]

For a given \(\chi\), if \(g\left(\chi, t\right) < 0\) then the contribution of state \(\chi\) is relatively negligible with respect to the maximally significant one, thus \(\chi\) can be discarded.

The interior of the Lieb-Robinson lightcone constists of the relativly significant states \(\chi\), and is defined by the equation

\[ g\left(\chi, t\right) > 0, \]

in analogy with the metric of special relativily.

Test of the apriori lightcone estimate on example problem.

Let us apply our lightcone estimate to the problem (5).

We just propagate \(g\left(\chi, t\right)\) forward in time and compute it for the chain sites. We add the latter when the metric becomes positive:

r_cut = 10**(-4)

ti = [0]
rho_plus = np.zeros((ns, ns), dtype = np.cdouble)
j = 1
for i in range(nt-1):
    # update rho_plus
    psi = tools.as_column_vector(psi_lc[:, i])
    rho_plus += tools.dyad(psi, psi) * dt
    # find max eigenvalue
    pi_max, _ = tools.find_largest_eigs(rho_plus, 1)
    # check whether we are inside the lightcone
    site_sig = rho_plus[j + 1, j + 1]
    #print(site_sig)
    #print(pi_max)
    #print('-------------')
    lr_metric = site_sig - r_cut * pi_max
    if lr_metric > 0:
        ti.append(i)
        j += 1
        if j==ns-1:
            # if all sites are coupled
            break

if ti[-1] < nt-1:
    ti.append(nt-1)
        

Let us visualize the resulting \(m\left(t\right)\):

       
x = []
y = []

for i in range(len(ti)-1):
    x.append(t[ti[i]])
    x.append(t[ti[i+1]])
    y.append(i+1)
    y.append(i+1)

plt.xlabel('time')
plt.ylabel('number of coupled modes')

plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7fd91a6f3160>]

On a smaller scale to see the step-wise growth:

plt.xlabel('time')
plt.ylabel('number of coupled modes')

plt.plot(x[:20], y[:20])

[<matplotlib.lines.Line2D at 0x7fd91a756b90>]

Let us solve the quench problem to test the apriori estimate. First, with all the sites coupled:

import secondquant as sq
hs_atom = sq.fock_space(num_modes = 1, max_total_occupation = 1, statistics = 'Bose')
m = 9 
n = 7 
fs_chain = sq.fock_space(num_modes = m, max_total_occupation = n, statistics = 'Bose') 
hs_joint = sq.fock_space_kron(hs_atom, fs_chain)
b_hat = hs_joint.annihilate
b_hat_dag = hs_joint.create
sigma_m = b_hat[0]
sigma_p = b_hat_dag[0]
sigma_x = hs_joint.sigmax(0)
a_hat = b_hat[1:]
a_hat_dag = b_hat_dag[1:]

Eat = 1.0
h = 0.05
e = 1.0

def f(t):
    return(0.1*np.cos(t))

Hconst = Eat * sigma_p @ sigma_m + h * sigma_p @ a_hat[0] + h * sigma_m @ a_hat_dag[0] \
    + e * sum([a_hat_dag[i] @ a_hat[i] for i in range(m)]) \
    + h * sum([a_hat_dag[i + 1] @ a_hat[i] + a_hat_dag[i] @ a_hat[i + 1] for i in range(m-1)])

def H(t):
    return Hconst + f(t) * sigma_x

def Hi(ti):
    return H(ti * dt + dt/2)

K = hs_joint.dimension
psi_ini = np.zeros(K, dtype = complex)
psi_ini[0] = 1 

tmax = 80
dt = 0.01
t = np.arange(start = 0, stop = tmax, step = dt)
nt = t.size

nq = np.zeros(nt)

for i, psi in tools.evolution(start_index = 0, end_index = nt, H = Hi, dt = dt, initial_state = psi_ini):
    nq[i] = (np.conj(psi) @ sigma_p @ sigma_m @ psi).real

Then let us solve the Schrodinger equation but now keeping only \(m\left(t\right)\) modes (from the apriori estimate):

nq2 = np.zeros(nt)

def eval():
    
    psi_begin = np.copy(psi_ini)
    
    for j in range(len(ti)-1):
        a = ti[j]
        b = ti[j+1]
    
        if a>nt-1:
            return
        b = min(b, nt-1)
    
        Hmconst = Eat * sigma_p @ sigma_m + h * sigma_p @ a_hat[0] + h * sigma_m @ a_hat_dag[0] \
            + e * sum([a_hat_dag[i] @ a_hat[i] for i in range(j)]) \
            + h * sum([a_hat_dag[i + 1] @ a_hat[i] + a_hat_dag[i] @ a_hat[i + 1] for i in range(j - 1)])
    
        def Hm(t):
            return Hmconst + f(t) * sigma_x
    
        def Hmi(ti):
            return H(ti * dt + dt/2)
    
        #nonlocal psi_begin
        for i, psi in tools.evolution(start_index = a, end_index = b, H = Hmi, dt = dt, initial_state = psi_begin):
            nq2[i] = (np.conj(psi) @ sigma_p @ sigma_m @ psi).real
            
            
        
        psi_begin = psi
        
eval()
              
plt.plot(t[:-2], nq2[:-2], label = 'only inside apriori lightcone')
plt.plot(t, nq, label = 'all sites')
plt.legend()
plt.xlabel('time')
plt.ylabel('occupation of qubit')

Text(0, 0.5, 'occupation of qubit')

We see that the apriory lightcone works!